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ABSTRACT  : 

This  paper  presents  recent  advances  in  a  computing  method  for  unsteady  inviscid  incompressible  flows 
around  thin  lifting  surfaces  with  arbitrary  motion.  This  involves  the  velocity  field,  deriving  from  a  poten¬ 
tial,  the  pressure  field  and  the  free  vortex  sheet  characteristics  (geometry  and  vortex  strength).  When  con¬ 
sidering  high  angles  of  attack,  the  separation  lines  along  which  the  vortex  sheets  are  shed  are  the  trailing 
edge,  as  well  as  wing  tips  or  parts  of  the  leading  edge.  Compared  to  classical  methods  of  singularity,  such 
as  panel  methods,  an  improved  vortex  wake  model  was  constructed.  It  consists  in  representing  the  wake  by 
a  vortex  particle  concept,  particularly  well  suited  for  the  prediction  of  its  unsteady  deformations.  This 
approach  stems  from  the  unsteady  vortex  sheet  theory  established  by  MJJDRY,  in  which  the  vortex  sheet  is 
considered  as  a  median  layer,  on  which  is  founded  the  concept  of  con  tin  uous  vortex  particle.  This  ensures 
the  computing  method  to  compare  favorably  with  the  classical  methods  in  terms  of  flexibility  and  com¬ 
puting  costs.  Computed  unsteady  flow-field  simulations  are  presented  for  several  cases  and  compared 
with  some  available  experimental  data. 


1  INTRODUCTION 

In  order  to  develop  computational  methods  to  analyse  complex  aircraft  flight  conditions,  the  simulation  of 
unsteady  aerodynamics  and  the  resulting  wake  dynamics  is  necessary.  It  is  much  needed  in  aeronautical 
disciplines  such  as  flight  dynamics  and  flight  simulation  and  for  time-dependent  structural  load  analysis. 
Detailed  solution  of  the  complete  nonlinear  fluid  dynamic  equation  along  time-dependent  flight  path  is 
still  costly,  particularly  because  of  the  computational  grid  to  encompass  large  wake  histories  along  com¬ 
plex  arbitrary  path.  An  alternative  approach  is  to  use  simplified  fluid  dynamic  equations  while  retaining 
the  three-dimensional  nature  of  both  a  wing  geometry  and  its  flight  path.  In  the  general  framework  of 
inviscid  incompressible  irrotational  fluid,  thin  wakes  are  represented  by  vortex  sheets.  Two  broad  types  of 
methods  currently  exist,  differing  in  the  way  they  consider  the  wake. 

The  first  considers  the  wake,  usually  approximated  by  trapezoidal  elements,  as  a  distribution  of 
doublets  on  its  surface.  Assuming  constant  doublet  density,  this  is  equivalent  to  closed  vortex  rings1  and 
corresponds  to  a  vortex  lattice,  which  moves  with  time  according  to  the  local  velocity  field.  Computer 
developments  have  allowed  direct  numerical  approach  of  such  three-dimensional  non-linearized  problems 
(Albano  and  Rodden2,  Djodjodihardjo  and  Widnall3,  Katz  and  Maskew4,  Mook5,  Johnson  and  al.6). 

The  second  type  of  methods  considers  the  wake  as  thick  enough  to  be  represented  by  volume  vor¬ 
tex  particles  (Rehbach7,  Morchoisne8,  Huberson9,  Winckelmans  and  Leonard10,  Buron11).  These  methods 
are  based  on  a  volume  discretization  of  the  vortex  vector,  i.e.  the  velocity  curl.  The  time  variation  of  the 
discrete  distribution  of  fluid  particles,  with  each  particle  supporting  this  vortex  vector,  is  described  with  a 
Lagrangian  integro-differential  formulation  based  on  the  Helmholtz  equation. 

Actually,  both  are  not  really  well  founded  in  theory  as  more  or  less  based  on  discrete  formulations. 
The  vortex  lattice  approach,  often  used  because  of  its  competitive  numerical  capabilities,  is  built  on  a  dis¬ 
crete  view  of  the  problem  and  overlooks  the  theoretical  problem  of  a  continuous  wake  evolution  with 
time,  so  that  this  can  lead  to  numerical  problems  for  large  wake  distortions.  The  above  vortex  particle 
models  allow  to  develop  an  accurate  deformation  scheme,  but  has  the  drawback  of  high-cost  numerical 
algorithms,  especially  for  the  Helmholtz  equation,  using  complex  regularized  schemes10'11. 
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This  work  presents  a  general  nonlinear  theory  for  time-dependent  problems  of  lifting  and/or  pro¬ 
pulsive  three-dimensional  airfoils  of  variable  geometry  in  arbitrary  motion  in  a  potential  flow.  In  many 
cases,  and  this  will  be  done  in  this  paper,  they  can  be  considered  as  thin  enough  to  be  assumed  to  be  thin 
surfaces.  The  airfoil  geometry  may  be  arbitrary  in  the  sense  that  leading,  tip  and  trailing  edges  may  be 
curved  or  kinked  and  the  airfoil  may  have  arbitrary  camber  and  twist. 

The  problem  was  first  formulated  in  theory12'13  with  no  reference  to  any  discrete  approach.  It 
mostly  relates  to  the  unsteady  thin  wake  model  established  by  Mudry14  .  An  exact  continuous  scheme  was 
developed  from  this.  It  particularly  focuses  on  the  wake  shedding  and  time  deformation,  so  as  to  develop 
an  efficient  computational  tool  for  the  unsteady  airloads  and  wake  geometry.  The  method  has  been  im¬ 
plemented  for  computation  of  helicopter  rotor  wake15.  For  such  applications,  it  is  sufficient  to  consider 
free  vortex  sheets  emanating  only  from  the  trailing  edge  because  rotor  blades  have  a  high  aspect  ratio. 

This  paper  presents  new  developments  of  the  method  for  low  aspect  ratio  wings  with  evolutive 
vortex  sheets  emanating  from  leading,  tip  and  trailing  edges. 


2  DESCRIPTION  OF  THE  METHOD 


2.1  Background 

Consider  A,  a  thin  airfoil  in  motion  in  a 
surrounding  fluid,  and  3  the  shedding  line 
of  its  wake  £  (figure  1).  In  the  limiting  case 
of  an  inviscid  flow,  the  wake,  resulting 
from  the  union  of  the  two  distinct  boundary 
layers  issuing  from  the  lower  and  upper 
surfaces  of  the  wing,  is  an  infinitely  thin 
layer.  It  is  a  slip  surface,  and,  as  it  cannot 
sustain  any  pressure  difference,  is  a  dis¬ 
continuity  surface  for  the  tangential  veloc¬ 
ity  component  of  the  fluid  particles  located 
on  each  side  of  this  surface. 


Fig.  1  -  Wake  in  inviscid  fluid 


The  fluid  is  inviscid  and  incompressible  over  the  entire  irrotational  flow  field,  excluding  wing  and  wake. 
Therefore,  a  velocity  potential  (|)  is  defined  in  an  inertial  frame  of  reference,  which  satisfies  Laplace’s 
equation,  and  is  submitted  to  the  following  boundary  conditions: 
flow  tangency  (zero  normal  velocity  across  the  wing), 

quiescent  free  stream  fluid  (automatically  fulfilled  with  the  below  problem  formulation), 
smooth  flow  off  the  shedding  edges  (Kutta-Joukowski  condition), 
pressure  continuity  across  the  wake. 

Green’s  theorem  allows  to  write  the  velocity  potential,  on  and  outside  the  surfaces  A  and  £  : 

(1)  PeIR3  ^(P)  =  -^-||o(M)^^<lS-^||(i(M)^^dS 

A  1  E 

where  o  and  p  are  doublet  strengths  over  the  surfaces,  respectively  A  and  £.  P  is  a  field  point,  nM  the 


— ^  ^ 

outward  normal  in  M,  and  r  =  MP ,  r  =  |f| .  For  P  on  respectively  A  or  £,  the  first  or  second  integral  in 

Equation(l)  is  a  Cauchy  singular  integral,  which  has  to  be  evaluated  in  the  principal  value  sense. 

The  derivation  of  the  velocity  potential  on  any  surface  S  at  any  point  in  the  field,  leads  to  the  resulting 
formulation  (Mudry14)  : 


(2) 


grad't>='ib/i’ 


(nM  AVp(M))Ar 


dS  -| - , 

4nA 


r:  /  h(M)gradM 


as 


AdM 


in  which  dS  is  the  oriented  edge  of  S,  V  is  the  surface  gradient  of  a  scalar  quantity  defined  on  S  rela¬ 


— > 


tively  to  a  local  frame  and  gradM  ,  the  gradient  of  a  scalar  quantity  defined  in  M  and  evaluated  at  a  point 
M  of  S.  This  expression  shows  the  usual  two  terms:  a  surface  integral,  involving  the  surface  gradient  of 
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the  doublet  density  distribution  on  the  surface,  and  a  curvilinear  one  involving  the  value  of  doublet  density 
distribution  along  the  edge  of  the  surface.  Mudry11  showed  that  the  second  term  in  (2)  has  to  be  equal  to 
zero,  so  that  p  has  to  be  zero  at  any  free  edge  of  S,  union  of  the  thin  airfoil  A  and  its  wake  E. 

2.2  Thin  wakes  and  the  concept  of  vortex  particles 

In  the  general  framework  of  potential  flows,  Mudry 16  developed  a  general  theory  for  unsteady  wakes,  con¬ 
sidered  as  vortex  sheets  (discontinuity  surface  for  the  tangential  velocity)  represented  by  classic  double 
layers.  Mudry  founded  his  theory  on  the  median  layer  concept17,  which  was  first  introduced  by  Helmholtz 
in  terms  of  median  velocity.  The  vortex  sheet  E  representing  the  wake  is  characterized  by  the  non-unique 
pair  of  functions: 

-  <15 ,  the  "median  parametrisation",  shown  in  figure  2, 
which  determines  the  geometrical  shape  of  the  vortex 
sheet  and  its  deformation  : 


(3)1  (q’q  )en<=IR  t  e  I  dR  x  =  S(q1,q2,t)e  Et 
x  =  (x1,x2,x3)eIR3 


Fig.  2  -  Parametrisation  of  a  surface. 


-  y  ,  the  associated  "median  vortex  density"  that  describes  the  vortex  sheet  strength.  It  is  a  function  of  the 
local  velocity  jump  and  of  the  geometrical  description  of  the  wake  carried  out  by  the  median  parametrisa¬ 
tion: 

f  ^  ->  \ 

dG5  3g5 

dq1  3q2 

,  J 

The  wake  convection  is  described  by  the  following  evolution  equation  that  G5  verifies  : 

3®(q‘-q2,t)  U+  +  U” 


(4) 


y  (q'-q2’  t)  =  na  [u]  = 


r-A  [dS  2  ds 

Alul=1V+1V 


■  =  U-(<3(q‘,q2,t)) 


(5) 

dt  2 

where  U  is  the  median  velocity  field  due  to  the  flow,  i.e.  the  half-sum  of  the  two  fluid  surfaces  velocities. 
The  position  of  a  vortex  particle  being  defined  by  the  point  N  :  1 1— >  ON  =  £D(q1,q2,t  j,  (G5,y) 

clearly  states  the  concept  of  continuous  vortex  particle  to  represent  the  wake. 

From  the  choice  of  the  parametrisation  0  together  with  the  pressure  continuity  condition  through 
the  wake,  Mudry  demonstrated  the  fundamental  time-conservation  property  for  the  two  contravariant 
components  f  (a=l,2)  of  the  median  vortex  density.  The  wake  shedding  is  taken  into  account  by  speci¬ 
fying  that  q1  is  the  shedding  edge  parameter  linked  to  a  shedding  edge  representation:  q,  and  q2  the  shed¬ 
ding  time:  x.  It  can  be  deduced  that  the  two  contravariant  components  of  the  vortex  density  are  determined 
once  for  all  when  the  particle  is  shed  with  the  determined  shedding  relative  velocity  at  and  relatively  to  3  : 
^30(q,x,x)j  3h(q,x) 


(6) 


V„  = 


at 


ye 


at 


ah . 


where  —  is  the  absolute  velocity  of  the  shedding  edge  points  by  using  a  lagrangian  representation 

at 

— ^  ^ 

OP  =  h  (q,t)  of  the  shedding  edge  3  in  the  inertial  frame. 

The  velocity  induced  by  the  vortex  sheet  E  takes  the  forms  : 


(7) 


a  ||x-0(q,x,t)||  47t3L  '  ||x-ffi(q,x,t)|| 

where  is  the  velocity  potential  and  3  means  shedding  edge.  The  integration  domain  is  a  known  range 
defined  by  the  definition  plane  (q,x)  of  the  parametrisation  (5  .  It  does  not  explicitly  depend  on  time. 
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The  problem  of  computing  the  wake,  which  consists  in  determining  the  position  and  density  of  the  vortex 
particles  once  shed,  reduces  to  solving  the  integro-differential  equation  that  governs  the  new  position  of 
the  vortex  particles  at  every  moment. 

2.3  General  Nonlinear  Equations 

The  unsteady  motion  of  the  airfoil,  pre¬ 
scribed  in  an  inertial  frame  of  reference 
Rg  =  (O  ;  X,  Y,  Z),  is  written  in  an  airfoil- 
fixed  frame  R  =  (O’  ;  x,  y,  z),  as  shown  in 
figure  3,  and  corresponds  to  a  rigid  airfoil 
motion  (velocity:  vE).  An  optional  relative 
motion  within  R  describes  the  deforma¬ 
tions  of  the  airfoil,  in  addition  to  its  aver¬ 
age  motion. 

Whereas  the  geometrical  shape  of  the  wake  Z  is  unknown,  it  is  known  for  the  thin  airfoil  A.  As  for 
the  wake,  a  parametrisation  %  defined  in  the  inertial  frame.  The  flow  tangency  condition  governs  the  po¬ 
tential  discontinuity  c|)+-  o’  =  [o]  on  A.  [d>]  has  to  be  equal  to  zero  on  the  borders  of  A  and  Z,  excluding 

their  common  border  3.  Through  3,  its  variations  will  be  governed  by  the  Kutta-Joukowsky  condition18. 

Figure  4  shows  how  the  physical  problem  is  then  formulated  in  two  domains,  respectively,  related 
to  the  definition  plane  of  the  parametrisations  G5  and  x  ■  This  is  one  of  the  original  features  of  this  ap¬ 
proach,  as  the  problem  is  formulated  and  solved  for  both  vortex  sheets  (airfoil  and  wake)  in  their  respec¬ 
tive  two-dimensional  parametrisation  planes. 


Fig.  3  -  Airfoil-fixed  frame  R  and  galilean  frame  Rc. 


O'M  =  x(e',e2,t)  O’N  =  05(q,T:,t) 

O'P  =  X  =  ro(q,T,t)  =  h(q,t) 


Fig.  4  -  Physical  field  and  parametrisation  planes. 
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The  nonlinear  theoretical  problem  of  a  thin  airfoil  in  an  unsteady  flow  is  settled.  Flow  tangency, 
shedding  relative  velocity,  Kutta-Joukowsky  and  wake  time-variation  equations  form  a  nonlinear  integro- 
differential  system  to  be  solved  which  solution  requires  a  numerical  approach. 

In  many  usual  configurations,  the  wake  shedding  at  a  trailing  edge  can  satisfactorily  be  prescribed 
on  the  shedding  edge  flight  path.  In  this  case,  Vre  is  assumed  to  be  equal  to  the  flight  path  velocity  of  3 


- ,  as  it  is  predominant.  Thus,  the  integro-differential  system  solution  is  simplified  since  Vre  is  no 

(It 

longer  an  unknown.  On  the  other  hand,  when  considering  high  angles  of  attack  and  thus  separated  flows 
with  wake  shedding  from  airfoil  tips  or  leading  edge,  this  assumption  is  no  longer  valid.  Indeed,  in  this 
case,  the  shedding  relative  velocity  Vre  has  to  consider  both  terms,  the  resulting  velocity  induced  by  the 


airfoil  and  its  wake 


dt 


<9h 

,  as  well  as  the  kinematic  velocity - .  In  these  conditions  Vre  is  an  unknown 

dt 


for  the  shedding  process,  which  will  have  to  be  included  in  an  iterative  computing  process. 


The  pressure  jump  can  be  computed  by  the  Bernoulli  equation  in  the  airfoil-fixed  system  of  refer¬ 
ence.  The  pressure  difference  across  the  airfoil  is  then: 

(8)  ^  =  (U*  +  VE)VM+3M 

p  v  '  dt 

The  pressure  and  the  potential  field,  allow  to  compute  the  normal,  lift  and  drag  forces  and  corresponding 
moments,  as  well  as  all  the  desired  flow  characteristics,  such  as  surface  velocity  surveys. 


2.4  Discrete  formulation  and  numerical  aspects 

Discretisation  of  the  two  geometrical  surfaces,  A  and  £,  is  performed  through  approximations  of  x  and 

G5  defined  in  relation  with  a  subdivision  of  their  respective  definition  planes.  The  simplest  possibility  is  to 
approximate  the  definition  planes  by  the  union  of  small  uniform  quadrilateral  panels  and  to  approximate 
X  and  G5  by  their  values  at  the  center  of  these  panels.  These  approximate  values  of  x  and  03  are  then  the 
control  or  collocation  points,  where  discrete  equations  are  solved. 

A  step  distribution  of  the  doublets  over  the  panels  has  been  retained  for  the  airfoil.  This  distribu¬ 
tion  is  equivalent  to  a  vortex  ring  distribution  in  the  definition  planes  of  x  ■  Consequently,  this  type  of 
approximated  distribution  is  not  very  different  from  a  classical  vortex  lattice  approach,  except  that  the 
discretisation  is  performed  in  the  parametrisation  planes. 

The  availability  of  the  continuous  theory  allowed  Mudry  to  introduce  for  the  wake  a  discrete  vor¬ 
tex  particle  concept  deriving  directly  from  the  its  continuous  formulation,  by  considering  a  step  distribu¬ 
tion  of  the  contravariant  components  fL  (cx=  1 ,2)  of  the  vortex  density  over  the  panels.  Using  i,j  as  indexing 
letters  referring  to  a  wake  panel,  the  velocity  induced  by  a  discrete  vortex  particle  distribution  is  then  : 


where  h1  and  h2  are  increments  of  the  discrete  plane  (q,x).  In  equation  9,  the  03-dcrivatcs  referred  to  the 
definition  plane  have  to  be  approximated,  using  a  finite-difference  scheme. 

It  is  well  known  that  any  discretisation  of  this  type  (vortex  ring,  discrete  vortex  particle)  intro¬ 
duces  a  fluid  velocity  singularity,  whereas  the  continuous  formulation  does  not.  This  is  the  case  as  x 
tends  towards  the  vortex  particle  point,  i.e.  when  r  tends  towards  zero.  In  order  to  avoid  numerical  diver¬ 
gence  and  thereby  regularize  the  singular  behavior  of  the  discrete  velocity  field,  a  regularization  process  is 
required.  According  to  reference  12,  what  has  been  retained  is  multiplying  the  induced  velocity  by  the 
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r 

coefficient  — - - ,  where  rc  is  known  as  the  regularization  radius.  In  many  study  tests,  this  radius  has 

r~+rc” 

been  chosen  equal  to  one  hundredth  of  the  mean  chord. 

Numerical  Method 

The  classical  numerical  scheme  is  based  on  a  time-stepping  procedure.  The  time  interval  [0,T]  is 
divided  into  intervals  At.  At  initial  instant  (t  =  0)  no  wake  exists.  At  each  time  step  the  airfoil  is  moved 
along  its  flight  path  and  the  corresponding  shed  wake  can  be  predicted  from  the  set  of  discrete  equations. 
At  each  time  step,  the  wake  strip  shed  during  the  current  time  step  is  a  row  of  free  vortex  ring  elements. 
Their  positions  have  to  be  calculated  by  the  shedding  process,  which  is  governed  by  two  coupled  equa¬ 
tions,  the  Kutta-Joukowsky  condition  and  the  time  variation  equation.  Each  panel  is  associated  with  a 
vortex  ring  or  a  vortex  particle  by  changing  the  discretisation  model. 

At  each  control  point  on  the  airfoil,  the  local  velocity  induced  by  the  airfoil  and  its  wake,  has  to 
satisfy  the  flow  tangency  condition.  At  this  stage,  the  unknowns  are  the  strengths  of  the  airfoil  vortex  rings 
and  the  characteristics  of  the  vortex  rings  currently  shed.  The  position  of  each  of  the  other  wake  discrete 
vortex  panicles  has  been  determined  by  solving  the  wake  time  valuation  equation  during  the  previous  time 
step,  and  its  strength  by  the  Kutta-Joukowsky  condition,  once  for  all  when  shed. 

In  the  case  of  a  trailing  edge  wake  shedding,  the  positions  of  these  vortex  rings  elements  are  de¬ 
rived  directly  from  the  shedding  assumption  Vre  =  -3h/3t .  Thus,  in  the  physical  domain,  the  shed  vortex 
ring  elements  correspond  to  the  interval  covered  by  the  trailing  edge  during  the  time  step.  Concerning 
their  strengths,  they  are  determined  by  applying  the  Kutta-Joukowsky  condition.  In  this  case,  the  discrete 
Kutta-Joukowsky  condition18  then  consists  in  setting  the  strengths  of  the  vortex  ring  elements  equal  on 
both  sides  (airfoil  and  wake)  of  the  trailing  edge.  Using  the  wake  shedding  procedure,  the  discretisation  of 
the  flow  tangency  equation  derived  at  the  control  points  yields  to  a  linear  system  of  algebraic  equations, 
where  the  unknowns  are  the  strengths  of  the  airfoil  vortex  rings,  and  which  take  into  account  the  strengths 
of  the  currently  shed  vortex  rings  because  of  the  Kutta-Joukowski  condition.  Further  details  on  the  proce¬ 
dure  may  be  found  in  reference  12  and  13.  At  each  time  step,  the  solution  is  obtained  by  inverting  the  ma¬ 
trix  of  the  linear  system.  If  the  shape  of  the  airfoil  remains  unchanged,  then  the  matrix  inversion  occurs 
only  once. 

When  considering  airfoil  tips  or  leading  edge  shedding,  each  vortex  ring  element  of  the  wake  ship 
is  shed  from  the  separation  line  according  to  a  local  shedding  relative  velocity  (Eq.  6).  This  velocity  de¬ 
pends  on  the  strengths  of  the  vortex  ring  elements  currently  being  shed,  which  are  unknown,  so  an  itera¬ 
tive  numerical  process  has  to  be  used.  At  each  iteration,  the  shedding  relative  velocity  is  obtained  from  the 
strengths  of  the  vortex  rings  in  process  of  shedding  and  computed  at  the  previous  iteration.  It  is  computed 
on  the  control  point  of  the  neighboring  airfoil  panels,  where  a  slip  condition  is  verified,  so  that  each  vortex 
element  is  shed  in  a  local  plane  tangent  to  the  airfoil  at  3 .  The  Kutta-Joukowsky  condition  and  the  flow 
tangency  equation,  which  have  now  to  be  applied  at  each  iteration,  are  processed  the  same  way  as  above. 

The  time  variation  equation  is  discretized  on  the  same  basis  and  its  numerical  solution  actually 
gives  the  new  position  of  the  vortex  particles  at  each  time  step.  This  solution  is  obtained  by  using  a 
Runge-Kutta  algorithm.  The  contravariant  components  of  y  are  computed  from  the  approximation  of  [eft] , 

which  has  been  determined  once  for  all  at  the  shedding  time  by  the  Kutta-Joukowsky  condition.  The 
resolution  order  of  the  RK-algorithm  depends  both  on  computation  time  and  accuracy  to  predict  the  wake 
geometry.  The  higher  the  order,  the  better  the  wake  rollup  is  simulated. 

One  part  of  this  work  was  to  carry  out  some  numerical  experiments  on  simple  configurations  in  order 
to  test  the  convergence  of  the  numerical  process  in  relation  with  the  discretisation  parameters,  such  as  the 
time-step  and  the  number  of  airfoil  panels  chordwise  and  spanwise.  The  results  have  shown  that  the  dif¬ 
ferent  numerical  schemes  behave  quite  well.  In  every  case,  convergence  was  obtained  by  improving  the 
discretisation  refinement  (increasing  the  number  of  airfoil  panels  chordwise  and  spanwise  and  decreasing 
the  time-step)  as  shown  in  references  12  and  13. 
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3  RESULTS 

3.1  Application  to  Flapping  Variable-geometry  Wings 

To  illustrate  the  capabilities  of  the  method,  we  present  the  calculation  of  flapping  wings  of  variable  ge¬ 
ometry,  based  on  works  about  the  forward  flight  of  birds16'17.  Consider  a  pair  of  flapping  wings,  which 
planform  shape  is  a  good  approximation  of  a  pigeon’s,  flying  with  a  constant  forward  velocity  U0  and 
hinged  about  a  longitudinal  axis  x  (figure  5). 


One  complete  period  of  the  flapping  motion  consists  of  a  downstroke,  an  upstroke  and  two  tran¬ 
sient  motions.  The  upstroke  and  the  downstroke  are  chosen  equal.  According  to  references  15  and  16,  the 
downstroke  produces  lift  and  thrust,  whereas  lift  and  drag  are  at  their  smallest  during  the  upstroke.  The 
assumption  of  constant  angular  velocity  of  flapping  (f)  is  applied  during  the  upstroke  and  the  downstroke. 
During  the  transient  motions,  (j>  has  a  sinusoidal  variation  .  A  wing  twist  allows  to  prevent  unrealistic 
geometrical  incidences  resulting  from  the  combination  of  the  forward  and  flapping  motions.  This  twist  is 
such  as  the  resulting  geometrical  incidence  does  not  exceed  a  set  incidence  qm  during  the  downstroke, 
and  is  zero  during  the  upstroke.  Numerical  tests  showed  that  the  best  results  were  obtained  with  a  rectilin¬ 
ear  twist  rotation  axis  set  at  the  three-quarter  line  chord.  The  following  results  were  obtained  with  the  fol¬ 
lowing  parameter  values  : 

-  10c 

reduced  flapping  period  :  T= - — ,  0  =  0°,  (f)max  =  80° ,  amax  =  20° 

U0 


Due  to  its  symmetry,  only  half  of  the  problem  is  considered.  Discretisation  for  the  wake  geometry 
computation,  involved  4  equally  spaced  chordwise  elements  and  30  spanwise  with  a  non-dimensional 

time-step  =  0.5 .  The  airload  computation  requires  a  more  refined  one  :  10  x  40  elements,  with 


— 5 —  =  0.2  .  For  2.5  periods,  and  without  any  particular  optimization,  the  computation  run-time  values 

Co 

range  respectively  from  some  minutes  to  a  few  hours  on  a  personal  computer.  The  numerical  results  show 
the  desired  characteristics,  i.e.  lift  and  thrust  are  essentially  produced  on  the  downstroke,  whereas  the  lift 
and  drag  are  negligible  on  the  upstroke  as  shown  in  figure  6. 


Fig.  6  -  Lift  and  drag/thrust  coefficient  time  variation. 
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Fig.  7  -  Wake  geometry  at 


^=20. 


3.2  Low  aspect  ratio  rectangular  wing 

Here  are  presented  results  of  computation  including  tip  wake  shedding.  It  concerns  the  flow  around  a  low 
aspect  ratio  (AR)  flat  plate  at  moderate  angles  of  attack. 

As  a  first  test  case,  the  sudden  acceleration  of  a  wing  that  was  initially  at  rest  is  investigated. 
Computed  results  for  lift  coefficient  and  spanwise  lift  distribution  are  shown  in  figures  8  and  9.  They  are 
compared  with  experimental  steady  measurements  obtained  by  Scholz  in  reference  21  for  an  AR  =  2  wing. 

U  t 

Values  are  obtained  for  a  reduced  time  of  — —  =  5  ,  which  corresponds  to  five-covered  mean  chord  c. 

c 


CL 


Fig.  8  -  Lift  coefficient. 


y/(b/2) 

Fig.  9  -  Spanwise  lift  distribution. 


The  wing  is  divided  into  10x20  elements,  chordwise  and  spanwise.  The  time  step  is  — - —  =  0.1 .  Com- 

c 

puted  results  with  the  current  method  compares  well  with  experimental  values.  In  order  to  predict  loads 
for  higher  angles  of  attack,  it  would  be  necessary  to  implement  a  leading-edge-emanating  wake.  Such 
model  is  being  implemented. 

The  second  test  case  presented  in  figure  10  concerns  a  unit  aspect  ratio  wing.  The  wing  is  divided 

U  At 

into  15x15  elements,  chordwise  and  spanwise.  The  time  step  is  — - —  =  0.06.  The  unsteady  motion  is  de- 

c 

fined  by  sinusoidal  oscillations  of  the  freestream  direction  with  time,  from  incidence  0°  to  ±14.5°.  Non- 
dimensional  period  of  oscillations  corresponds  to  one  covered  chord.  Incidence,  normal  load  coefficient 
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CN  and  “steady  “  normal  load  CNS  coefficient  are  represented  during  two  periods  on  figure  8.  CNS  is 
equal  to  CN  without  the  unsteady  contribution  of  the  term  —  . 


Fig.  10  -  Normal  force  coefficient.  Sinusoidal  oscillation  of  incident  freestream  (±14.5°). 


It  can  be  noted  that  the  unsteady  effect  of  forces,  proportional  to  — ,  leads  to  a  lead  phase.  CN  higher 

3t 

values  are  obtained  before  incidence  corresponding  values.  These  computed  results  compares  well  with 
results  obtained  by  Rehbach  in  a  similar  test  case  exposed  in  reference  7. 


A  pitching  oscillation  motion  of  the 
wing  is  also  computed.  Pitch  oscilla¬ 
tions  from  0°  to  ±11.5°  are  about  the 
leading  edge  and  the  non-dimensional 
period  of  oscillations  corresponds  to 
2.5-covered  chord.  The  difference  with 
the  above  motion  is  in  the  fact  that 
pitching  angular  velocity,  0,  is  taken 
account  in  the  motion  description  in 
the  wing-fixed  frame  defined  in  sec¬ 
tion  2.3.  It  contributes  to  emphasize 
unsteady  effects  on  loads  and  wake 
geometry.  Figure  1 1  shows  time  varia¬ 
tions  of  pitch,  lift  coefficient  and  pitch 
moment  coefficient. 


Fig.  1 1  -  Pitching  oscillations:  lift  and  pitch  moment  coefficient. 


Figure  12  displays  the  wake  geometry  obtained  at  the  reduced  time  step  —  =  3.5 .  Desired  tip  roll-ups 

c 

are  clearly  noted  and  stability  of  wake  simulation  is  demonstrated. 
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Fig.  12  -  Wake  geometry  at 


c 


3.3  Delta  wings 

Since  most  fighter  aircraft  incorporate  a  delta-wing  surface,  there  is  great  interest  in  studying  unsteady 
delta-wing  manoeuvers  including  pitching  to  large  incidence  angle  and  high-rate  rolling  motion.  Simula¬ 
tion  of  flow  separation  from  delta  wings  is  being  implemented.  In  order  to  predict  the  high-speed  flow 
induced  by  leading-edge  vortices,  the  current  method  in  its  state,  and  especially  for  delta  wings,  needs  to 
include  a  more  refined  shedding  process,  which  is  being  testing. 

First  results  are  presented  in  figure  13  for  an  unit  aspect  ratio  delta  wing  with  an  unit  mean  chord. 
This  is  obtained  for  an  impulsive  start  at  25  degrees.  Figure  13  displays  the  wake  geometry  obtained  at  the 

reduced  time  step  =  4.  The  reduced  time  step  is  ^°°At  =  0.04  and  the  wing  model  consists  of  80  pan- 
c  c 

els  with  8  chordwise  equally  spaced  panels. 


iu 


=  4. 


Fig  13  -  Wake  geometry  at  t 
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4  CONCLUSIONS 

On  a  basis  of  a  rigorous  and  complete  continuous  theoretical  formulation,  a  numerical  method  is  devel¬ 
oped  to  solve  the  unsteady  3D  aerodynamics  for  thin  lifting  system.  Applications  of  the  current  method 
exposed  in  this  paper  demonstrates  that  the  theoretical  model  exposed  in  section  2  can  extend  to  separation 
lines,  along  which  the  vortex  sheets  arc  shed,  such  as  wing  tips  or  parts  of  the  leading  edge. 

The  theory  leads  to  a  view  of  the  vortex  wake  similar  to  that  of  a  classical  vortex  lattice  approach, 
but  uses  a  discrete  vortex  particle  concept,  which  is  particularly  well  suited  for  computation  of  the  un¬ 
steady  deforming  wake.  In  fact,  the  present  method  allows  a  treatment  of  wake  deformations  comparable 
with  the  vortex  particle  methods  quoted  as  references,  though  without  their  limitations,  and  compares  with 
classical  vortex  lattice  approaches  in  terms  of  computing  costs.  For  all  results  presented  in  this  paper,  the 
computation  run-time  doesn’t  exceed  one  hour  on  a  computer  with  a  Pentium  II  350Mhz  processor. 

The  method  in  its  current  state  is  based  on  the  assumption  that  the  flow  is  symmetrical.  In  order  to 
simulate  more  complex  and  no  symmetrical  flight  configurations,  first  future  developments  need  to  free 
this  assumption  and  also  to  consider  multi-element  systems. 
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PAPER:  10 

AUTHOR:  DR.  LEROY 

Question  by  Dr.  Hummel:  Have  you  applied  your  method  to  a  flow  situation  for  which  it  is  known  from 
experiment  that  vortex  breakdown  occurs? 

Answer:  Not  all  flow  situations  are  covered. 

Question  by  Dr.  Nangia:  About  numerical  stability  -  the  close  approach  problem  of  2  vortex 
lines  or  more  -  Do  you  assume  cut-off  distances  to  stabilize  the  solutions? 

Answer:  We  assume  cut-off  distances,  which  directly  derive  from  the  discretization  models  of 
such  methods.  For  the  current  method,  we  explained,  in  our  corresponding  paper,  that  these 
distances  have  been  chosen  to  be  equal  to  one  hundredth  of  the  mean  chord  in  many  study  tests. 

Question  2:  Fidelity  of  vortex  flows  and  roll-ups  -  in  60's  and  70's  this  method  was  pioneered. 

The  main  difficulty  was  that  rolls-ups  in  close  detail  were  not  accurate  enough,  e.g.  for  delta 
wings.  This  is  presumably  due  to  discretization  distances.  How  fine  does  one  need  to  go? 

Answer:  Because  of  computer  capabilities,  we  have  now  the  possibility  of  improving  the 
discretization  refinement  (increasing  the  number  of  airfoil  panels  chordwise  and  spanwise  and 
decreasing  the  time-step)  in  order  to  obtain  "convergence"  of  the  discretization  model  in  relation 
with  the  discretization  parameters.  Some  examples  are  shown  in  our  paper  and  some  quoted 
among  the  references.  But  in  complex  configurations,  e.g.  delta  wings,  in  order  to  predict  the 
high-speed  flow  induced  by  leading-edge  vortices,  it  is  also  necessary  to  implement  the  fine 
shedding  process,  particularly  for  vortex  sheets  emanating  from  leading  edges,  and  to  increase  the 
resolution  order  of  the  Runge-Kutta  algorithm  as  well.  The  higher  the  order,  the  better  the  wake 
rollup  is  simulated.  This  is  what  we  are  testing  for  delta  wing  computations  with  the  current 
method. 

Question  by  Dr.  Lamar:  Are  there  comparisons  of  this  method  with  other  methods,  including 
that  of  Professor  Eddie.  Lan,  University  of  Kansas,  published  in  1981? 

Answer:  We  present  some  comparisons  with  other  methods  (vortex  lattice  and  vortex  particle 
classical  methods)  in  our  paper  and  in  references  12  and  13,  but  we  have  no  comparison  with 
Professor  Lan's  method  because  we  do  not  know  yet  of  this  reference. 
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